Can the Ising critical behavior survive in 
'\b'[ non-equilibrium synchronous cellular 

O ■ automata? 

(N 



t^ ! Kazumasa Takeuchi * 

^ . Department of Physics, The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo, 



B 

C/5 



CJ 



X 



113-0033, Japan 



Abstract 



C^ , Universality classes of Ising-like phase transitions are investigated in series of two- 

dimensional synchronously updated probabilistic cellular automata (PC As), whose 
,._i ■ time evolution rules are either of Glauber type or of majority-vote type, and degrees 

r^ . of anisotropy are varied. Although early works showed that coupled map lattices 

O \ and PCAs with synchronously updating rules belong to a universality class distinct 

from the Ising class, careful calculations reveal that synchronous Glauber PCAs 
should be categorized into the Ising class, regardless of the degree of anisotropy. 
^ \ Majority-vote PCAs for the system size considered yield exponents v which are 

0^ ' between those of the two classes, closer to the Ising value, with slight dependence 

t_^ . on the anisotropy. The results indicate that the Ising critical behavior is robust with 

■^ i respect to anisotropy and synchronism for those types of non-equilibrium PCAs. 

CD ' There are no longer any PCAs known to belong to the non-Ising class. 

O 
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1 Introduction 



The concept of universality class, which has been undoubtedly one of the cen- 
tral issues in equilibrium physics, is widely believed to hold in some range 
of non-equilibrium systems [1,2]. That is, even far from equilibrium, many 
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microscopic details of systems become irrelevant at transition points, and a 
set of critical exponents depends only on a small number of basic, macro- 
scopic ingredients. The "basic ingredients" comprise, e.g., spatial dimension, 
symmetries and conservation laws, as in equilibrium systems, but then it is 
natural to ask "do there exist any additional relevant parameters intrinsic to 
non-equilibrium?" 



The Ising universality class is also observed in non-equilibrium systems. For 
example, probabilistic cellular automata (PCAs) and coupled map lattices 
(CMLs) with up-down symmetry can exhibit Ising-like transitions by varying 
a control parameter such as a coupling constant. Grinstein et al. [3] showed 
that, suppose a coarse-grained dynamics of such systems is described by a 
Langevin equation, irreversibility due to the broken detailed balance is irrele- 
vant and the considered models should fall in the Ising class, using the standard 
dynamic renormalization group treatment. This prediction has been actually 
confirmed in various kinetic Ising models with asynchronously updating rules 
[4]. However, Marcq et al. [5] numerically found that Ising-like transitions of 
some two-dimensional non-equilibrium CMLs separate into two distinct uni- 
versality classes: the Ising class, and a new universality class, called "non-Ising 
class" hereafter, where the correlation length exponent v is 0.90(2) [here num- 
ber (s) between parentheses indicate the uncertainty in the last digit (s) of the 
quantity], which differs from the Ising value z/ = 1. Ratios of exponents l3/v 
and 7/z/ are common to the two classes, namely 0.125 and 1.75, respectively. 
Since CMLs with synchronously updating rules form the new class, whereas 
asynchronously updated ones belong to the Ising class, synchronism is thought 
to be a relevant parameter. After their work, several synchronous systems such 
as stochastic CMLs [6] , a logistic CML at the onset of a non-trivial collective 
behavior [7], and even PCAs [8,9] have been investigated and the existence of 
the non-Ising class has been observed again, aside from a few exceptions [6,10]. 
Thus the importance of synchronism, which may be related to the existence 
of an external clock, has been attracted much attention [11]. 



However, it is obvious that the synchronous updating does not immediately 
bring about the non-Ising critical behavior: there exist synchronous PCAs 
which respect the detailed balance and therefore we can safely say that they 
belong to the Ising class, thanks to the equilibrium universality hypothesis. 
For example, an isotropic PCA with a Glauber transition rate satisfies that 
condition. On the other hand, a Glauber PCA with completely anisotropic 
interaction was numerically studied by Makowiec and Gnacihski (MG) [9], and 
concluded to be in the non-Ising class. The above two observations naturally 
lead to a supposition that a degree of anisotropy may affect the selection of 
the universality classes. The aim of this paper is, therefore, to examine the 
relation between anisotropy and the universality classes. 



2 Models 



Two series of PCAs with up-down symmetry, namely Glauber PCAs and 
majority- vote PCAs, are numerically investigated. Both of them are on a two- 
dimensional, square lattice of size L x L with periodic boundary conditions. 
A local variable, or "spin," s*j G { — 1,+1} is assigned to each lattice point, 
where indices i and j denote Cartesian coordinates, and t is the discrete time. 
Each site {i,j) is simultaneously updated according to a specific local tran- 
sition probability Pij{±l\{s}) that the spin takes a value ±1 after one time 
step from a spin configuration {s}. For the Glauber PCAs, it is defined as 

Pij{±l\{s}) = -{1 ± tanhg[sij + Si+ij + Sij^i + a{si-ij + Si,j+i)]}, (1) 

where g is a coupling constant acting as a control paramter, and < a < 1 
indicates the degree of anisotropy. We can realize a variety of Glauber PCAs 
with different degrees of anisotropy by varying the value of a. An isotropic 
case corresponds to a = 1, in which it is easily shown that the detailed bal- 
ance is satisfied with stationary distribution P({s}) = ^ e~^^^^^\T-C{{s}) = 
— J2i,j In cosh (7 {sij + Sj-ij + Sj+ij + Sjj-i + Sjj+i) , where H({s}) is an ef- 
fective Hamiltonian and 1/Z is a normalization constant. Since the Hamil- 
tonian respects the up-down symmetry and consists only of short-range in- 
teractions, the equilibrium universality hypothesis asserts that this isotropic 
Glauber PGA should fall in the Ising class despite the synchronous updating 
scheme. On the other hand, the detailed balance does not hold for a 7^ 1. In 
particular, the completely anisotropic case a = is already reported by MG 
to have the exponent u = 0.93(3) and to belong to the non-lsing class [9]. 

The other series of PCAs is that of majority-vote PCAs, defined as 

Pij{±l\{s}) = -{1 ±gsgn[sij + Si+ij + Sij_i + a(si-ij + Sij+i)]}- (2) 

For this model, the degree of anisotropy is classified into only 3 types: com- 
pletely anisotropic < a < 0.5, intermediate a = 0.5, and isotropic 0.5 < 
a < 1. The detailed balance is violated in all of them [3,12]. The completely 
anisotropic case is well-known as the Toom PGA [13], whose exponent was 
also estimated by MG to be z/ = 0.86(4), i.e. the non-lsing value [8,?]. 



3 Methodology 



All of the simulations are implemented on lattices of size up to L = 96 by the 
following procedure. We start from random initial conditions and discard first 
to = 10^ time steps as transients. It is sufficiently long to consider the systems 



in the time-asymptotic attractor, since the correlation time is Xcorr ~ 5 x 10^ 
in the worst case. Then time series of the averaged spin m\ = {1/ L"^) J2i,j sjj 
is used to calculate the magnetization M^, higher order moments M|" , the 
susceptibility xl and the Binder's cumulant Ul [14], defined by 



M^ = (K|), Ml") = (|mir), 

XL = L\M^^-Ml), f/^(^) = 1 - -3-i^, (3) 






where (■ ■ ■) denotes the expectation value obtained by integrating an observ- 
able during the duration T, such as (|m^|) = (1/T) X)iiio+i l^^il- Absolute 
values of m^ are taken in eqs. (3) as usual, since finite-size effects allow sign 
reversals of the instantaneous magnetization even in the ordered phase. The 
integration time T is 8 x 10^ for the isotropic Glauber PGA and 1.5 x 10^ for 
the others, and thus much longer than both the correlation time Tcorr ^ 5 x 10^ 
and the sign reversal time r^ev ^ 5 x 10^. Therefore, the quantities in eqs. (3) 
are capable of representing the corresponding ensemble averages with good 
accuracy. 



4 Measurement of critical exponents 



Gritical exponents /?, 7, v are estimated for the Glauber PGAs with a = 0, 
0.25, 0.5, 0.75, 1 and the majority-vote PGAs with a = 0, 0.5, 1. The way 
of measurement is almost similar to that of Marcq et al. and later works 
[5,6,8,9,10], that is, to exploit finite-size scaling laws in equilibrium which 
are empirically known to hold in the non-equilibrium Ising-like transitions. 
The following shows the process of measurement for the Glauber PGA with 
a = 0.75 as a typical case. 

The first to do is to find the critical point g = g^- In order to locate it, we 
adopt the standard method using the Binder's cumulant UL{g) [14]. Since 
the cumulant has the scahng form of Uiig) = U{{g — gcjL^^'^), it becomes 
independent of L at criticality, i.e. Uiigc) = U* for all L. The quantity U* is 
also a universal number. Figure 1 shows plots of UL{g) and their polynomial 
fitting curves for various system sizes. 

Focussing on intersections between them, a drift toward larger values of g is 
clearly observed (the inset of fig. 1). This is caused by remaining contribu- 
tions of irrelevant operators, which exist in finite-size systems [14]. Taking it 
into consideration, the scaling function of the cumulant should be modified to 
Uiig) = U{{g — gc)L^^'^, OL~^)^ where only one irrelevant term with a relax- 
ation exponent uj is assumed to be dominant. An expansion of the modified 
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Fig. 1. Plots of the Binder's cumu- 
lant UL{g) for the Glauber PGA with 
a = 0.75. System sizes are, from the 
smallest slope, L = 16, 20, 24, 28, 32, 
40, 48, 64, 80, 96. Symbols correspond to 
raw data and curves indicate 5th order 
polynomial fits. The inset is a magnifica- 
tion of the intersection region. 



Fig. 2. Plot of coordinates g of intersec- 
tions between curves ULi{g) and 1/12(9) ■ 
The intersections are located by fitting 
raw data with polynomials of orders 4 
to 9, and superimposed on the same fig- 
ure. The range of the orders is chosen so 
that the fitting curves adequately trace 
the raw data. The broken line indicates 
the linear regression, whose y-intercept 
yields an estimate of Qc- See eq. (4). 



scaling form in powers of (g — gc)L^^'^ and L "^ up to the 1st order yields the 
coordinate of the intersection between Ul^^q) and 1/12(9), namely, 



g = gc + A 
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where A and B are expansion coefficients. Consequently, by plotting g with 
respect to (L^"^ — L2'^)/(Li — L2) for fixed, reasonable values of uj and z/, 
and by searching for uj which gives the highest correlation coefficient for it, we 
obtain estimates of both uj and gc- Figure 2 is a plot for u = 1 and the best 
value of UJ, where we can clearly see the linear dependence, which gives the 
coefficient of determination r^ = 0.909. The critical point gc is then estimated 
as a y-intercept of the regression line, gc = 0.34981_|3J in this case. Here, the 
superscript (subscript) indicates the confidence interval in the plus (minus) 
direction, which is evaluated as the region with sufficiently large coefficient 
of determination, namely, (1 — r^) < (1 — max^^ r^) x 1.1. The estimation of 
U* is carried out quite similarly, which results in U* = 0.6107_^3J. This is 
in good agreement with the 2D-Ising value U* = 0.611(1) [15]. A problem of 
this method may be that it requires the value of u, which is measured later, 
but the difference in the estimates with respect to the assumed value of z/, 
1.0 or 0.9, turns out to be quite subtle: 2 x 10"^ for gc and 2 x 10"^ for U*, 
i.e. negligible. The mentioned improvement is therefore useful and expected 
to lead to more reliable estimation of critical exponents. 



Now we proceed to the measurement of critical exponents. The following finite- 



size scaling relations are made use of to achieve it: 

a,lnMiU~LV^ 9,lnMf|,^~LV^ 

9,lnMf)|,, ~Ll/^ dgUL\g^^L^'\ 

V, ^ [m]y[m'] ~ Ll/^ (5a) 

and 

M^(^,) ~ L-'^/^ mP(^,)V2^^-/3/.^ 

m1^)(^,)i/^ ~ L-^/^ XL(^e) ~ i^^/^ (5b) 

where [m"] = <9gM)" |g^, and the derivatives are evaluated by using polyno- 
mial fittings of appropriate orders. Irrelevant operators can affect the scaling 
relations (5) again. We, therefore, perform the measurement as follows. (I) 
First, we plot a quantity in eqs. (5) in the log-log scale, using a polynomial 
fit of a fixed order, and check the linear dependence. (Ila) If the contributions 
of irrelevant fields are already suppressed for the smallest size considered, 
L = 16, a simple linear fit giving the lowest chi-square in the log- log plot 
yields an estimate. (lib) Otherwise, finite-size corrections are employed to 
achieve an asymptotic value of the slope, after Marcq et al. and subsequent 
works [5,6]. An expansion is made similar to that used for U^lg) above, e.g. 
L^^^'^dglnMLlg^ ~ Ci + C2L~^ 1 followed by searching for the best oj and v 
(or other exponents), in the sense of lowest chi-square. This method, how- 
ever, does not work sometimes due to rapid convergence of correction terms, 
or statistical errors in the raw data. In that case, (lie) we use a linear fit in 
the log-log plot neglecting data points subject to the finite-size effect. Care is 
taken to ensure that the observed scaling behavior does indeed correspond to 
the asymptotic regime, in all of the cases. Finally, (III) we repeat the aforesaid 
procedure for polynomial fits of different orders in (I), and also for all of the 
quantities in eqs. (5) which give the same exponent. 

A value of v for the Glauber PC A with a = 0.75 is thereby evaluated to 
be z/ = 0.988_[j^5N, as is shown in fig. 3. The main sources of errors are the 
uncertainty in the estimate of the critical point Qc and statistical errors due to 
the finite-time sampling. The former is estimated from the confidence interval 
of Qc, while the latter is from the dependence of exponents on the quantity 
used for the finite-size scaling, and on the degree of the fitting polynomials. 
Note that the final range of errors is given as the union of each error region 
considered, i.e. as wide as possible. Recalling the value of z/ for the Ising 
and non-Ising class, 1 and 0.90(2) respectively, we can clearly conclude that 
the Glauber PGA with a = 0.75 belongs to the Ising class. Values of jS/u 
and 7/// are obtained in the same way, which are P/u = 0.1271_/28) and 
7/z/ = 1.749_L4J, and again completely consistent with the Ising values. 
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Fig. 3. (a) Scaling relations of the 7 quantities in eq. (5a) for the Glauber PCA with 



a = 0.75, where X 
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from top to bottom. Corresponding slopes are 1.017, 1.013, 1.011, 1.014, 1.012, 
1.011 and 1.006, respectively, (b) Deviations from the scaling laws, which are found 
to be faint. Note that the scale in the y-axis is identical for all of the subplots. 

5 Results and discussions 



Critical exponents of all of the models considered are measured in the same 
manner, except for the following two points, (i) Systems of size L = 16, 20, 24, 
28, 32, 40, 48, 64, 80, 96 are examined for the Glauber PCAs with a = 0.25 
and 0.75, while those of size L = 16, 24, 32, 40, 48, 64, 80, 96 are considered 
for the others, since finite-size corrections are appropriately caught by them, 
(ii) Since no finite-size effects are observed in intersections of U^lg) for the 
Glauber PCA with a = 0, values and errors of Qc and U* are determined by 
means and standard deviations of the coordinates of the intersections. 



The measured exponents are summarized in table 1. Our results for the Glauber 
PCAs clearly show that they fall into the Ising class. While the isotropic 
case a = 1, which is at equilibrium, should actually be there, the results 
for all the anisotropic Glauber PCAs are quite unexpected since they are 
non-equilibrium, synchronous PCAs. In particular, we obtain a value of the 



exponent for the completely anisotropic case a = as u - 
which conflicts with the estimate by MG, z/ = 0.93(3) [9]. 



0.992 



+(28) 
-(26) 



(fig- 4), 



On the other hand, the results for the majority-vote PCAs are rather unclear: 
the estimated exponents u are between the Ising and non-lsing value. They 
are closer to the Ising exponent and this tendency slightly increases with 
isotropy. There seem to be two possible interpretations for it. One is that the 
critical exponents of the majority-vote PCAs depend continuously on system 
parameters, anisotropy in our demonstrations. This is, however, commonly 
understood as an exceptional case, at least in equilibrium. Instead, it is rather 
natural to consider that what we see here is a part of an extremely slow 



Table 1 

Critical points and exponents of the PCAs examined. 



model a 


gc U* P/u 7/1/ u 


1 

0.75 
Glauber 0.5 
0.25 



0.31173+g; 0.6098+jJJ) 0.1271+JJ2 l-757+|5Jj 0.993+[Jg 
0.34981+g; 0.6107+gJ 0.1271+gj l-749+[;;j 0.988+[;^J 
0.40407+!S 0.6100+g°) 0.1254+;2! 1-757+!!) 0.994+!?J! 

-(3) -(7) -(42) -(26) -(15) 

0.49049+[3j 0.6097+[2J O.I266+J2JJ l-753+[JJj 0.994+gj 
0.65855+[Jgj 0.6123+jJgj O.I239+J32} l-753+[j2J 0.992+J2gj 


1 
majority- vote 0.5 



0.73151+!^! 0.6105+S 0.1247+!i^| 1-753+^q^) 0.979+^^^^ 

-(4) -(4) -(17) -(9) -(17) 

0.783324t|iJ) 0.6104+J5J O.I26O+JJ2) l-753+|7J 0.962+JJqJ 
0.82248+g5 0.6140+JJJj O.I249+I24J l-740+j;°j 0.956+[Jg 


Ising [15] 
non-Ising [5,6] 


0.611(1) 0.125 1.75 1 

w 0.611 « 0.125 «1.75 0.90(2) 



convergence in finite-size scalings toward the Ising asymptotic behavior, which 
is estimated to be reached for L ~ O(IO^). The shght dependence on the degree 
of anisotropy is then related to the difference in relaxation constants, which 
may be attributed to coupling intensity between sites, and/or an additional 
length scale caused by some hidden coherent structures such as in [16], if they 
exist. In any case, further studies are essential to give a conclusion on it. In 
addition, the estimate for the completely anisotropic majority-vote PCA, or 
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Fig. 4. Measurement of the critical ex- 
ponent V for the Glauber PCA with 
a = 0. Each symbol and line corresponds 
to X = a^lnM^^I^,, 9plnMi')|g„ 
OglnMilg^, T/2, dgUilg,, V4, Ve from top 
to bottom. Slopes are 1.023, 1.016, 1.012, 
1.005, 1.005, 1.005 and 1.005, respec- 
tively. 



Fig. 5. Measurement of the critical expo- 
nent u for the majority-vote PCA with 
a = 0. Each symbol and line corresponds 
to X = dg\nM^^^\g^, SglnMf Ig,, 
dglnML\g„ V2, V4, Ve, dgUL\g, from top 
to bottom. Slopes are 1.046, 1.044, 1.041, 
1.044, 1.041, 1.039 and 1.026, respec- 
tively. 



the Toom PCA, v = 0.956_L J (fig. 5) is again incompatible with the value 
from MG, u = 0.86(4) [8,9]. 

The discrepancies between our estimates and MG's are caused by a few differ- 
ences in steps toward estimation. One is in the sampling method. We started 
from random initial conditions and sampled T = 1.5 x 10^ consecutive data 
after discarding first to = 10^ time steps, whereas MG chose ordered initial 
conditions, in which all spins are +1, and set to ^ 10^, T = 10^ as the price for 
repeating independent simulations not more than 5500 times. Since to and T 
of MG are in many cases comparable with the correlation time Tcorr ^ 5 x 10^, 
and comparable with or much shorter than the sign reversal time Trcv ^ 5 x 10^, 
we consider that the sampling in MG is statistically insufficient and the in- 
fluence from the ordered initial conditions may remain. Another origin of the 
discrepancies is a way to determine transition points Qc- We estimated them 
solely from the crossings of the cumulant Uiig), with finite-size corrections 
if possible, while MG first made a guess from the crossings and then deter- 
mined them by searching for g^ which gives the ratios of exponents /9/z/ and 
7/z/ identical to those of the Ising class. This is, however, rather perilous since 
estimates of critical exponents are very sensitive to various errors, including 
statistical ones. Moreover, in general, a small uncertainty in critical points Qc 
leads to much larger errors in critical exponents. In fact, if we assume the 
value of Qc in MG for the Glauber PCA with a = 0, namely Qc = 0.6580, we 
reproduce their result u ~ 0.93. For the reasons above, we believe our results 
are more reliable. 

In conclusion, we have investigated the Ising-like phase transitions in syn- 
chronous PCAs, namely the Glauber PCAs and the majority-vote PCAs, 
with different degrees of anisotropy. Our calculations reveal that the Glauber 
PCAs belong to the Ising class regardless of the degree of anisotropy, and the 
majority-vote PCAs are also expected to do so, though the latter remains to 
be clearly shown. The results indicate that the Ising critical behavior is robust 
with respect to anisotropy and synchronism for the PCAs, which coincide with 
the theoretical prediction of Grinstein et al. [3] and observations by Sastre and 
Perez where some degree of deterministic dynamics is required to bring about 
the non-Ising critical behavior [6]. 
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